#### illustration of confounding ####

# create data
set.seed(12)
mydat <- data.frame(x=NULL,z=NULL,y=NULL)
for(i in 1:5) {
   x<-(2.5*(5-i)):(12+2.5*(5-i))
      y<-2*i+.25*x+rnorm(13,sd=.8)
      z<-rep(letters[i],13)
      temp.dat<-data.frame(x,z,y)
      mydat<-rbind(mydat,temp.dat)
  }

plot(y~x,data=mydat, cex=.8)
abline(lm(y~x,data=mydat), col='grey70')

# an important variable has been ignored
points(mydat$x,mydat$y,col=as.numeric(mydat$z), pch=14+as.numeric(mydat$z))

# separate regressions by level of z
for(i in 1:5){
   out <- lm(y~x,data=mydat[mydat$z==letters[i],])
   lines(mydat[mydat$z==letters[i],"x"], predict(out),col=i)
}

out1<-lm(y~x+z,data=mydat)
for(i in 1:5){
   lines(mydat[mydat$z==letters[i],"x"], predict(out1)[mydat$z==letters[i]] , lty=2)
}
legend('bottomleft',c('parallel','separate','ignoring groups'), lty=c(2,1,1), col=c(1,1, 'grey70'), bty='n',cex=.8)

## Table 1: interaction ###
anova(lm(y~x*z,data=mydat))

# Table 2: ignoring covariate
summary(lm(y~x,data=mydat))
# confounding
summary(lm(y~x+z,data=mydat))